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Abstract 

Micro-plasticity theories and models are suitable to explain and predict mechanical response of devices on length scales 
where the influence of the carrier of plastic deformation - the dislocations - cannot be neglected or completely averaged 
out. To consider these effects without resolving each single dislocation a large variety of continuum descriptions has 
been developed, amongst which the higher-dimensional continuum dislocation dynamics (hdCDD) theory by Hochrainer 
et al. (Phil. Mag. 87, pp. 1261-1282) takes a different, statistical approach and contains information that are usually 
only contained in discrete dislocation models. We present a concise formulation of hdCDD in a general single-crystal 
plasticity context together with a discontinuous Galerkin scheme for the numerical implementation which we evaluate 
by numerical examples; a thin film under tensile and shear loads. We study the influence of different realistic boundary 
conditions and demonstrate that dislocation fluxes and their lines’ curvature are key features in small-scale plasticity. 
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1. Introduction 


Plastic deformation of metals has been utilized by man since the copper age and the knowledge 
of how to process metals (by e.g. forging) has constantly grown. However, the physical mechanisms 
underlying the empirical procedures could not be understood until the crystalline structure of metals 
was investigated by Rutherford at the beginning of the 20th century. Subsequent attempts to explain 
the discrepancy between the theoretically predicted shear strength of a metal and the experimentally 
observed yield stresses lead to the concept of the ‘dislocation’ - a linear crystal defect - which was 
proposed in the 1930s independently by Orowan (19341), Polanyi (1934) and Taylor (1934a). Two 


decades later Kondo (1952), Nye (1953), Bilby et al. ([1955 ) and Kroner (1958) independently in¬ 


troduced equivalent measures for the average plastic deformation state of a crystal in the form of a 
second-rank dislocation density tensor. This ’Kroner-Nye tensor’ is introduced to link the microscopi¬ 
cally discontinuous to a macroscopically continuous deformation state and is the fundamental quantity 
in Kroner’s continuum theory of dislocations. It has been used widely until today (see e.g. the re¬ 
lated works in Acharya and Fressengeas (2012); Taupin et al. (2013); Le and Gunther ( |2014 )). This 
tensor, however, only captures inhomogeneous plastic deformation states associated with so-called 
geometrically necessary dislocations (GNDs) and does not account for the accumulation of so-called 
statistically stored dislocations (SSDs) in homogeneous plasticity. This limits the applicability of the 
classical dislocation density measure within continuum theories of plasticity. 
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Phenomenological continuum models for plasticity which are not based on dislocation mechanics 
have been successful in a wide range of engineering applications. They operate on length scales where 
the properties of materials and systems are scale invariant. The scale-invariance, however, breaks 
down at dimensions below a few micro-meters, which is also a scale of growing technological interest. 
These microstructural effects become more and more pronounced in small systems and lead to so- 
called ’size effects’ (e.g. Ashby, [1970 


Arzt, 1998; Stolken and Evans, 1998; Greer and De Hosson 


2011). Phenomenological continuum theories incorporate internal length scales by introducing strain 
gradient terms - sometimes based on the consideration of GND densities - into their constitutive 


Fleck et ah. 

1994; 

Nix and Gao, 

1998; 

Gurtin, 

2002; 

Gao and Huang, 

2003; 

Zhang 


et al., 2014|) but are not able to consider fluxes of dislocations or the conversion of SSDs into GNDs 


and vice versa. Refined continuum formulations (with or without strain gradients) include additional 


information in a mechanism-based approach (Engels et ah, 

2012; 

Li et aLf |2014|) 

or take multi-scale 

approaches by directly including information from lower scale models (Wallin et al., 

2008; 

Xiong et al., 


2012 ). 


Discrete dislocation dynamics (DDD) models (e.g. Kubin and Ganova, 1992 


1997; Fivel et ah, 19971 


et al., 2007; Zhou et al. 


Devincre and Kubin 


Ghoniem et al.[[20001 [Weygand et al. , 2002; Bulatov and Gai , 2002; Arsenlis 


2010 


Po et ah, 2014) contain very detailed information about the dislocation 


microstructure and the interaction and evolution of dislocations and have been successful over the last 
two decades in predicting plasticity at the micro-meter scale. DDD simnlations allow to investigate 
complex plastic deformation mechanisms but are, however, due to their high computational cost 
limited to small system sizes/small densities. 

A different approach which is closely related to DDD and which generalizes the classical continuum 
theory of dislocations was undertaken by Groma et al. (Groma, 1997; Groma et ah, 2003). They used 
methods from statistical physics to describe systems of positive and negative straight edge disloca¬ 
tions in analogy to densities of charged point particles. The resulting evolution equations are able to 
faithfully describe fluxes of signed edge dislocations and the conversion of SSDs into GNDs (and vice 
versa). This approach has been successfully used by a number of groups (e.g. [Yefimov et al.l |20C 


Kratochvil et ah, 

2007; 

Hirschberger et ah, 

2011; 

Scardia et ah. 

2014) 


A generalization to systems 
of curved dislocation loops, however, is not straightforward. Pioneering steps into that direction have 
been undertaken by Kosevich (1979); El-Azab (2000); Sedlacek et al. (2003). Furthermore, ’screw- 


edge’ representations have been introduced as an approximation by Arsenlis et al. (2004); Zaiser and 


Hochrainer (2006); Reuber et al. (2014) and Leung et al. (2015); Xiang and co-workers developed 


a model for the evolution of curved systems of geometrically necessary dislocations (Xiang, 2009). 


Their model also includes line tension effects and was e.g. applied to model Frank-Read sources (Zhu 


et ah, 2014). A new approach based on statistical averages of differential geometrical formulations 


of dislocation lines has been done by Hochrainer ( 

Hochrainer, 

2006; 

Hochrainer et ah. 

2007; 

Sandfeld 


et ah, 2010) who generalized the statistical approach of Groma towards systems of dislocations with 


arbitrary line orientation and line curvature introducing the higher-dimensional Continuum Disloca¬ 
tion Dynamics (hdGDD) theory. The key idea of hdGDD is based on mapping spatial, parameterized 
dislocation lines into a higher-dimensional configuration space, which contains the local line orienta¬ 
tion as additional information. This is particularly useful in situations with complex microstructure 
(Sandfeld et ah, 2010, 2015). In order to avoid the high computational cost of the higher-dimensional 
configuration space, ’integrated’ variants of hdGDD - denoted by CDD - have also been developed 
recently and their simplifying assumptions already have been benchmarked for a number of situation 


(Hochrainer et al. 

,2009; 

Sandfeld et al., 

2011; 

Hochrainer et al., 

2014; 

Zaiser and Sandfeld, 

2014; 

Mon- 

avari et ah. 

2014) 

. Furthermore, very recently one GDD variant has been coupled to a strain gradient 


plasticity model (Wulfinghoff and Bohlke, 2015) and was used to study size effects of a composite 
material. Until now hdGDD nonetheless serves as reference method for all GDD formulations, since 
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it can be considered as an almost exact continnnm representation of ensembles of curved dislocations 
and allows to access many relevant microstrnctnral information. 

In this article, we derive and formulate the governing equations for hdCDD in a general crystal 
plasticity context. We start by introducing the formulation for single-crystal plasticity and the elasto- 
plastic boundary value problem and connect this to the dislocation system in a staggered scheme. The 
evolntion of the dislocation density and the dislocation cnrvature density is compnted in representative 
slip planes depending on the dislocation velocity. This is approximated by a discontinnous Galerkin 
(DG) scheme introduced in Sect. with a local Fourier ansatz suitable for the higher dimensional 
configuration space. Numerical examples in Sect. |^- tension and shearing of a thin film - demonstrate 
the influence of passivated and non-passivated surfaces on the microstructure evolution. Furthermore 
we study the influence of the lines’ cnrvature on the plastic deformation behavior and also on the 
stress-strain response of the specimen and compare with the formnlation of Groma. 


2. Crystal elasto-plasticity based on the hdCDD theory 

In this section we introduce the elastic eigenstrain problem and its variational formulation (Sect.[^ 


and Sect. 2 . 2 ) and two different models for representing slip planes (Sect. 2.3) in a continnnm frame¬ 
work. Then we ontline the hdGDD theory extending the classical quantities from Krbner’s continnnm 
theory: the governing eqnations for the evolntion of dislocation density and curvature density as de¬ 
fined by the higher dimensional continuum dislocation dynamics theory are introduced in Sect. |2.5t our 
model for the dislocation velocity governing the dislocation dynamics is determined by the constitutive 
setting explained in Sect. 113 


2.1. A continuum model for single-crystal plasticity 

Let the reference configuration B he a bounded Lipschitz domain in and let d^B U d]<iB = dB 
be a non-overlapping decomposition into Dirichlet boundary OdB and Neumann boundary &^B. The 
position of a material point is denoted by x and the displacement of the body from its reference 
configuration by u(x, t). The distortion tensor 

Du = /3®^ + /3P^ (1) 

is decomposed additively into elastic and plastic parts and /3 p\ respectively. Small deformations 
are assumed so that the infinitesimal strain is given by 


e = £(u) = sym(Du). 


( 2 ) 


Plastic slip is assumed to take place on N slip systems determined by a unit normal and 
slip direction d* = Tb* on the s-th system, where b^ is the Burgers vector of length bg = |bs|. As 
an example, the face-centered cubic (FGG) crystal has A = 12 slip systems. In special situations 
symmetry can be exploited and the case A = 1 or A = 2 can be considered over the crystal of thin 
film via shearing and bending situations. 

The plastic shear strain in the slip system s is denoted by %. In the single crystal, we assume 
that the plastic part of the displacement gradient is given by the sum over the contributions from all 
active slip systems 


/3P' = (3) 

where M* = Ab^ 0 rUs = dg 0 rUs is the projection tensor accounting for the orientation of the slip 
system s. Depending on the vector of plastic shear strains 7 = ( 71 ,..., the plastic strain is 
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given by 


gPi = = sym(/3P') = y , Mf” = sym(M,). 


( 4 ) 


This defines the elastic strain 


S=> = £"(u. 7) = £(u) - £->(7) . (5) 

2.2. Variational balance equations 

The macroscopic equilibrium equation is given by 

— diver = fs in B, (6) 

with the body force fg and the constitutive relation for the Cauchy stress tensor 

o- = C : = C : (e - (7) 

depending on the elasticity tensor C. The macroscopic boundary conditions are 

u = Ud on OdB, ern = on d^B , (8) 


where ud is a prescribed displacement for Dirichlet boundary and gN is an applied traction for Neu¬ 
mann boundary. Let {7 = {u G E^): u = 0 on and assume that ud extends to B. Then, 

for given £p' = £^*( 7 ) , we have in weak form: find u G Ud + N such that 




• hu dx + 



gN • (5uda, 


SueU. 


( 9 ) 


This is complemented by an evolution equation for the plastic shear strain (the Orowan relation 
(Orowan, 1940)) depending on the dislocation density and the dislocation velocity. 


2.3. Dislocation densities, plastic shear strain, and the representation of physical slip planes 

Since continuum dislocation theories operate with averaged measures one has to consider slip planes 
in a consistent sense (i.e., consistent with the averaging). In the slip system s, we use a discrete set 
of ’crystallographic’ slip planes (SP) of distance As > 0 

Ts,g = { 7 .s,g + Tjds + X m3 : ( 77 , ^) G E^} , (10) 

where Zg^g G Tg^g denotes the origin of the local (77,^) coordinate system which is aligned such that 
the Burgers vector points into positive 77 direction and ^ points into the line direction of a positive 
edge dislocation; a position in the slip plane is denoted by r G Tg ^,. Each SP is expanded to a thin 
layer of width h < As (collecting a small number of physical slip planes) 

Bs,g = |zGB:z = r + with r G Tg^g and |C| < h/sj . (11) 

In our model, the dislocation density in the slip system s is represented by the average pg^g in the layer 
Bg^g such that the sum of all pg^g integrated over each layer equals the total dislocation fine length in B. 
Since the evolution of the dislocation density ps^g and the Orowan relation of the plastic shear strain 
are evaluated only in the crystallographic slip planes Ps^^, the continuum approach requires to extend 
the values to the body B. For this purpose, we introduce the orthogonal projection Pg^g-. B —Pg^j,, 
and for r G Pg^g \ B the plastic shear strain jg^g is extended by constant continuation. We consider 
two cases (compare Fig. [^: 
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Experimental observation of surface 
slip traces. This illustrative example 
shows a plastically deformed (cadmium) 
single crystal (height 800pm). 

Permission granted by DoITPoMS 



The physical model idealizes the real 
specimen: The ’crystallographic’ slip 

planes are evenly spaced with distance As 
which is typically much smaller than the 
simulated distance in the numerical model 
below. 



Numerical model 

The ’numerical slip planes’ represent a 
number of physical slip planes either in 
a quasi-discrete setting (Case 1, left) or in 
a fully continuous setting (Case 2, right). 



Case 1: direct SP representation Case 2: averaged SP representation 


Figure 1: From the real system to the numerical model: representation of crystallographic slip planes 
for one slip system. 


Case 1: Direct representation of crystallographic SPs. 

We set 


^7s,g(r) r G Bs,g for some g , 
0 else 


( 12 ) 


where the factor of As/h accounts for the fact that the plastic strain will be defined for the elastic 
BVP only in the blue regions in Case 1. The objective of Case 1 is to analyze how to choose As 
and h for our density-based micro-structure representation. This case also serves as a reference 
for Case 2 and was set up following the idea of discrete dislocation dynamics simulations: for 
the limit h —>■ 0 we can approach a discrete SP representation; this provides benchmark data for 
Case 2. 


Case 2: Averaged representation of crystallographic SPs. 

Alternatively, we average over multiple crystallographical SPs in order to arrive at a represen¬ 
tation of, e.g., dislocation density or plastic strain in which they are distributed field quantities 
- not only within the SP but also in direction of the slip system normal. Therefore, we first 
collapse a number of crystallographic SPs that are contained within a region of width As^h 
into one representative SP by summing up the respective dislocation field variables over As. The 
representative SPs are numbered by g. All points r of this layer belong to the domain 

Bs,g = |z e z = Tg-|-Cms with r-g E Pg^g and |C| < As/2|. (13) 

The domains Bg^g are non-overlapping with (J Bs;g = B. For the plastic shear strain in Bs,g we 
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define the average 


As ^> 

7 .,^ = — 2 ^ 


A5 


75,5f 


rs,gnBcBs^g 

At points between the representative slip planes 

A5-C C r n 

r ~ ^s,g^ + ~^Ps,g-\-l^ 5 C ^ [0? As] 

we define the plastic shear strain in the body B by linear interpolation 

7*(r) = ^^^7s,giPs,9r) + P7s,g+l{Ps,g+ir) . 


(14) 


(15) 


2.4- Kroner’s continuum dislocation dynamics theory 

The hdCDD theory extends Kroner’s continuum dislocation model (Kroner, 1958), which is based 
on the dislocation density tensor (the so-called Kroner-Nye tensor) 

a = V X /3P^ (16) 

relating the plastic distortion in an averaging volume to a dislocation densitjQ Thus, we have 

V-a = 0 (17) 


which implies that dislocation lines cannot end or start within the crystal. 

In the special case that all dislocations in an averaging volume form smooth bundles of non¬ 
intersecting lines which all have the same line orientation, they are geometrically necessary; the 
respective dislocation density is referred to as the GND density Based on this volume density, the 
dislocation density tensor depends on the average line orientation 1 and the Burgers vector b and is 
given by 


OL = G b . (18) 

Now we consider a specific slip system s. If all GNDs move with the same velocity = Usb x d^, one 
obtains the evolution equation 


dtOLs = -V X (vs X as). 


(19) 


The scalar velocity Vg = |vs| is a constitutive function depending on, e.g., stresses of the solution of 
the elastic boundary value problem 


2.5. A higher-dimensional model for continuum dislocation dynamics of curved dislocations 

In order to represent systems of curved dislocations and their evolution in a representative slip plane 
Hochrainer generalized Kroner’s theory (Kroner, 1958) and the statistical approach of Groma 


(Groma, 1997; Groma et ah, 2003) towards systems of dislocations with arbitrary line orientation 


and line curvature introducing the higher-dimensional Continuum Dislocation Dynamics (hdGDD) 
theory (Hochrainer, 2006; Hochrainer et ah, 2007; Sandfeld et ah, 2010). hdGDD is based on mapping 
spatial, parameterized dislocation lines into a higher-dimensional configuration space T^^p x S^, where 


^We note, that this density measure is dependent on the spatial resolution with which the plastic distortion is defined 


(Sandfeld et ah, 2010) 
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= M/ 27 rZ = [0, 27r) is the orientation space containing the local line orientation as additional 
information. The continuum representation of lines in this conhguration space requires the notion of 
a so-called generalized line direction L^, ^ and generalized velocity Ys,g together with the dislocation 
density tensor of second order cx^^g, which is also defined in the configuration space. This density tensor 
contains the Kroner-Nye tensor as a special case but is furthermore able to describe the evolution of 
very general systems of curved dislocations with arbitrary orientation; in particular the common 
differentiation between GND and SSD density becomes dispensable. Similar to Krbner’s or Groma’s 
frameworks, this continuum theory again also describes only the kinematics, i.e., the evolution of 
dislocation density in a given velocity field. The additionally available information of hdCDD, e.g. 
line orientation and curvature, however, is crucial for determining dislocation interaction stresses and 
modeling physically-based boundary conditions in a realistic manner. 

In the following, r = {r],^) is a point in the slip plane and (r, (/p) denotes a point in Tg ^, x S^. 
The dislocation density on T^^g must be understood as a volume density, and thus to obtain the total 
line length in the SP we have to integrate p^ g over Bs,g. Let !«((/?) = cosiy^dg -|- siniy^d^ x mg be 
the canonical spatial line direction, and Ls^g(r, cp) = (l^, ks^g)~^ defines the generalized line direction in 
the higher-order configuration space, where kg^g is the average line curvature. The dislocation density 
tensor of second order is given by 

V^)Ls,s(r, v?) ® b,, (20) 

and the evolution equation for this tensor has the form 

dtOilgir,(p) = - V X {\g^g{r,p) x a“g(r, 9 ?)) , (21) 


with V = (d^, d^, d^p), where the vector \s,g = {—Vs,gdpls, —Jjs,g'Yvs,g) denotes the generalized velocity 
in configuration space. For detailed information on derivations and implications of these equations 


refer to (Hochrainer et ah, 2007; Sandfeld et ah, 2010). 


Since the dislocation density tensor of second order is determined by the dislocation density and 
the average line curvature, inserting (20) in the evolution equation (21) results into the equivalent 

this takes the form 


system for Ps^g and kg^g. 


Introducing the curvature density Qg^g = Ps,gks, 


gi 


^tPs,g Y • (^Ps,gYg^g^ Qs,g^s,g ; (22a) 

P^tQs,g ^ ■ i.Qs,g'^s,g^ Ps,g(^s,g ' ^(Lg^g ' YVg^g^'^ (22b) 


(see (Ebrahimi et ah, 2014) for an interpretation of these two equations as integrated total line length 
in Bg^g and the number of closed dislocation loops, respectively). The system (22) is complemented 
by boundary conditions which define, e.g., whether a dislocation can leave the crystal (free surface) 
or not (impenetrable surface). In the latter case, the density flux through the boundary is set to zero. 


e.g. 


Ps,g'^s,g 


• n = 0, where n is the outward unit normal at the boundary dB. 


2.6. Averaged quantities in the slip planes 

From a formal point of view the Kroner-Nye framework with equations (18) and (19) corresponds 
to the higher-dimensional Hochrainer framework (20) and (21) where b and v* are replaced by their 
higher-dimensional counterparts. One can retrieve the Kr5ner-Nye tensor for the slip plane s from 
the density function by 





Ps,3(r,9?)b((/p) ® d(/p. 


(23) 
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Other classical measures can be derived as well, e.g. the total scalar density is obtained by integrating 
Ps,g{v^ip) over all orientations and the GND density is obtained as the norm of the GND vector 
Ks^g = [( 0 : 5 , 3 ) 11 , ( 05 , 3 ) 12 ]. The GND vector and the GND vector rotated by 7 r /2 derive as 

^27r ^27r 

>^s,g{r)= I ps,g{r,p>)\s{^)Aip, Kj^^(r) = / ps,g{r,p>)d^\s{‘p)A^, (24) 

JO Jo 


where dg,\s{^) = — sin ip d 5 +cos 9 ? dg x mg is the orthogonal line direction, and the total scalar densities 
are given by 

n27T ^27r 

ql°g{r)= qs,gir,p) dp. (25) 

Jo Jo 

The plastic shear strain is determined integrating the Orowan relation over the orientation space 


^t7.,3(r) 


p27r 

= bj ^ 

JO 




(26) 


Thus, for the case that the dislocation velocity is independent of the orientation this simplifies to the 
classical plastic strain rate equation dt'js,g = bsPl^gVg^g 


Orowan 


1940 


2.7. The dislocation velocity 

Continuum dislocation models are ’kinematic’ theories predicting the flux of density depending on 
the dislocation velocity Vg^g as a constitutive ingredient. Hence, in Kroner’s, Groma’s and Hochrainer’s 
theories alike, stresses from dislocation interactions are not a priori included in these theories and 
need to be determined separately. How to derive physically meaningful dislocation interaction stress 
components is a topic of ongoing research; details about rigorous analysis of some dislocation systems 


in a continuum framework can be found in (e.g. 

Zaiser et ah, 2001 

; Groma et ah, 

2003; 

El-Azab et ah, 

2007; 

Sandfeld et ah, 

2013; 

Schulz et al.| [2014; 

Scardia et ah, 2014 

). Here, we base the dynamics of 


our dislocation systems on the following assumptions for the velocity function: 


1. The scalar velocity Vg^g in r 5,3 is assumed to depend linearly on the stresses acting on dislocations. 

2 . The velocity function is decomposed into different stress contributions of two separate classes: 
the projection of the resolved stress Tg^g = M 5 : cr computed from the solution of the elastic 
BVP Q, and stresses governing short-range elastic dislocation interactions (the back stress 
line tension and the yield stress rig). 

3. We assume that dislocations move only if the yield stress is exceeded by the other stress contri¬ 
butions T^g = Tg^g — Tg g — rfg, and the flow direction is determined by the sign of 

With these assumptions the velocity function takes the form 


V 


s,9 


R sgn(r° )(|r° I-r]; ) if 






> 


0 


otherwise, 


(27) 


which assumes an overdamped dislocation motion neglecting inertial effects. H > 0 is the respective 
drag coefficient. 


^Note that it is the total density which governs the plastic strain rate and not the GND density; the latter would not 
be able to account for plastic strain due to expansion of a statistically homogeneous distribution of dislocation loops. 





















































The resolved stress Tgn depends on the global boundary value problem and includes the contribution 


from the eigenstrain, which governs the long-range interaction between dislocations (see e.g. Krbner 


1955; El-Azab et ah, 2007; Sandfeld et ah, 2013). 


The back stress is an approximation for the repelling forces between parallel dislocations with the 


same line direction (see e.g. 

Groma et al., 

2003; 

Hirschberger et ah, 

2011 ; 

Schulz et ak, 

2014 

). For a 

system of curved dislocations we adapt a formulation which was suggested and implemented in the 

related works ( 

Hochrainer 

,|2006; 

Sandfeld 

, 2010 

), and we define 




^tot 

Hs,g 


V • n 


s,g 


(28) 


This implies that the back stress at r in direction of ls( 9 ?) is proportional to the gradient of GNDs 
perpendicular to the line direction. Here, fi is the shear modulus, D e [0.4, 1 ] is a constant and was 


obtained in Groma et al. (2003) by statistical analysis of discrete dislocation systems. 

The line tension rib describes the self-interaction of a dislocation loop (a loop subjected to no 


s,g 


other stress would contract due to the line tension). In general this depends on the local character of 
the line (e.g., whether it is a screw or edge dislocation, cf. Foreman ( 1967| )). For simplicity we use an 
approximation by a constant line tension which is independent of the line orientation 


T** = 
s,g 


TsQ. 


,tot 


Hs,g 


(29) 


where the constant Tg G [Q.bjih'l, TO/ife^] describes the strength of the interaction and is the average 
curvature density. Finally, the yield stress is governed by a Taylor-type term in the case of only one 
active slip system 


_ 
s,g 


= a/ibg 


(30) 


with a constant a G [0.2,0.4] , see e.g. Taylor] (1934b) or the overview in Basinski and Basinski (|1979|). 


In general, the yield stress can depend on the dislocation density of all slip systems (Kubin et al. 


2008); this will be considered as well in Study 3 below. 


3. A numerical scheme for the reduced 2D system 

For the numerical evaluation of the physical behavior of our model we consider a 2D reduction of 
the fully coupled system, where we assume to have a homogeneous distribution over the direction 
and slip plane normals m^ in the x — y plane. This leads to drjPs,g = 5 ? 7 ?s ,3 = ^ 77^,3 = 0 everywhere in 
the system. In this section we describe the discretization for a single slip plane T = for simplicity 
we skip the indices s and g. 


Reduced hdCDD over ID slip plane. Rewriting (22) using drjP = drjq = = 0 yields 

dtp = —V • (pV) -I- qv 

dtq = —V • (gV) — p cos (pd(^ ^ cos ipd(^v -\- kd^v^ — q^^p ^ cos (pd(y -\- kd^v 


(31a) 

(31b) 


with V = and V = (usimp, — coscpd^u — kd^^vY. Assuming that the velocity does not 

exhibit any angular anisotropy the system (|M|) reduces further to 


dtP = -V • (pV) qv 

dtq = —V • (f?V) — p(cos (p)^5|u -|- q sin ipd^v . 


9 


(32a) 

(32b) 





















































(33) 


We rewrite (32) in the equivalent form of a linear conservation law for w = (p, g)^, i.e., 

dfW + V • F(w) + Bw = 0 in r X 
where V • F(w) = 9^(Fiw) + 9(^(F2w) with 


Fi = 


V sin ip 0 
0 V sin if 


p ^ f-cosipd^v 0 

^ 0 — cos ifd^v 


B 


0 —V 

(cos (pY&jv — sin (pd^v 


depending on the dislocation velocity v. 

For this first-order system we now construct a discontinuous Galerkin discretization in space, 
since this discretization is well adapted to conservation laws, leads to a better approximation for 
discontinuous weak solutions and is significantly less diffusive than standard continuous Lagrange 
elements (Hesthaven and Warburton, 2008). 


The Runge-Kutta discontinuous Galerkin (RKDG) method. The system (33) reads as follows: find 
w : r X 5”^ X [0, T] solving 


dtw + Aw = 0 in r X 


subject to the initial condition w(r, ip, 0) = Wo(r,(p), periodic boundary conditions w(^,p,0,t) = 
w(^, Tj, 2%, t) in ip, and where Aw = V • F(w) + Bw = Fi5^w + F29(^w + (5^Fi + c?^F 2 + B)w. 

In the first step we derive a semi-discrete discontinuous Galerkin scheme. Let Th = {K} be a 
triangulation of B, and assume that this triangulation is aligned with the slip plane r, i.e., r - U/e.Fr f 
with faces / G Jr = {dK fl F: K £ Th}- Let Zh be the set of all vertices of the triangulation. For a 
fixed face / = conv{zj_i, Zj} with Zj_i,Zj G F we observe 


A(w,V^)/xS'i = [ (V-F(w)-hBw) 

JfxS^ 


tp d^dcp 


/ ( — F(w) • Vi/?-|-Bw ■ r/?) d^dp-|- / n • F(w) ■ i/? 

JfxS^ Js^ 


Zj-l 


dp 


for all smooth functions ip: f x — > M^. For the discretization we choose an ansatz space Xf 
on every face defining a discontinuous ansatz space Xh = {iph G L 2 (F,M^): iph\f £ Xf} and the 
numerical flux 


= {{F/V’U}/,. + q^|[V>k]|/..s 

on the face intersections e = df H dfe (which are single points in a ID slip plane). Here, we choose 
the stabilization constant Cf = |V|, and the average and the jump along a normal n oriented from / 
to /e are given by 

{{FWh}}f,e ■-= ^ (Fw/,1/ + FWhlu'^ > [ W]/,e := W/,1/ - Wh\f, , 

respectively. On the open boundary, we set ^h\fe = ^h\f and w/j|^^ = —Wh\f for the impenetrable 
boundary. Now, defining locally 


{Af^h^h,'iph)fxs^ = 


/ {-F{wh) ■Xiph + 'Bwh-iph)(i^(iT+ n-F}^{wh) • iph 

//xSi Js^ 


dp 


Zj _ 1 
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yields the discrete operator by {AhMVh,'iJ)h)rxS^ = '^f{-^f,h^h,'^h)fxs^- We choose a DG ansatz 
space with Fourier basis functions Xf = Fk® F„, where = span{l, ..., are polynomials, and 
the truncated Fourier space is given by = span{l, cos((p), sin((y9),..., cos(n(y9), sin(n(p)} . Thus, the 
components of w^, = {ph,kh) have the form 

k n 

Phii, = XI + X cos{mip) + bim sin(m(/p))^ 

1=0 m=l 

and similar for qh. The Runge-Kutta time discretization is now obtained by the method of lines. 
Therefore, choosing a basis ..., of Xh yields the matrix formulation 


Mc}fw(t) + A w(t) = 0 


(34) 


with M = •0™'))^ ^ and A = This yields for the time step from tn to 

=w^- AtM”^A(w" - XtM"^A(w” - XtM"^A(w’^ - XtM'^A w”))) 

2 3 4 

for the classical explicit Runge-Kutta schemes of order 4 (see 
time integration methods in combination with DG schemes). 

Finite element discretization of the solid. Let 14 = {v G v\k G Pi(W) for K G Th} be a 

standard finite element space for the displacements and set 14(ud) = {v G 14^ v(z) = ud(z) for all 
nodal points z G ZhOdDBY For u G 14 the strain e(u) and the stress cr = C : £ is piecewise constant 
in K. Now, the coupled algorithm is defined as follows: 


Hochbruck et ah (2015) for alternative 


50) Select At > 0, tmax > 0, set n = 0, set initial values for on F^ ^. 

51) Set tn = nAt, ug = UD(tn), fe = hYn) and = gN(tn)- 
Compute 7 ” in B from 7 ”^ in F^^g (depending on Case 1 or 2 ). 

52) Evaluate the plastic strain £p'’"'|a' = Xs Ts and compute u” G 14(uo) with 

f £(u"') : C : £((5u) dx = f eP'’” : C : £((5u) dx + f fg-hudxT f gN-(^uda, hu G 14(0). 
Jb Jb Jb Jdt^B 

53) On f = K n K' C Fs,g set t”^|/ = and compute velocities vfg 

(eqn. (p7|)). 

54) Compute independently on every Fg^p by M explicit Runge-Kutta steps for 

dM] ) with step size At/M and fixed velocity vf g. 

55) If tn < tmax, Set u = u + 1 and go to SI). 


Since this scheme is in step S4) fully explicit, the Courant-Friedrich-Levy condition requires sufficiently 
small time steps. Also the coupling with the boundary value problem S2) is explicit. In our numerical 
tests we choose the global time step At and the local time step At/M small enough to observe 
convergence by comparing the results with different mesh resolution h and time steps At. 
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4. Numerical experiments 


We evaluate our model for a single-crystal thin film with idealized passivated and non-passivated 
surfaces in tensile and shear test settings. These are well established model tests for crystal plasticity 


models (see e.g. Liu et al. 

, 2011 

; Schwarz et ah, 2005| 

Fredriksson and Gudmundson 

, 2005 

; Fertig and Baker 


2009). The novelty of hdCDD, however, is 


that we can directly link the dislocation microstructure in almost DDD-like details to the macroscopic 
response. In the following, we study in particular the influence of the line curvature and two different 
physical boundary conditions in single- and multislip configurations. Additionally, we numerically 
evaluate effects due to averaging for Case 1 and 2. We note, that especially the line curvature is an 
important physical quantity that, however, only can be represented by few continuum models (e.g 

a fact that makes detailed comparisons with other approaches 


Sedlacek et ah, 

2003; 

Xiang, 

2009) 


difficult. 



(a) Tensile test: Displacements are prescribed at left and right (b) Shear test: Displacements are prescribed at bot- 
boundaries, top and bottom boundaries are free surfaces. tom and top boundaries, left and right boundaries 

are free surfaces. 


Figure 2: Geometry and boundary conditions of the investigated model systems for the tensile test 
(Study 1 and Study 2) and the shear test (Study 3). 


Geometry and slip system. We consider the configurations as shown in Fig. [ 2 ] for the investigation 
of the deformation behavior of a thin, single-crystalline A1 film assuming plane strain. The film is 
represented by a 2-dimensional body B = (0,10/) x (0,/) with / = 1 pm. We consider one or two 
active slip systems (A = 1 or 2) with 1-dimensional (crystallographic) slip planes determined by 

= cos 4> ei + sin 0 62 , = — sin (f) ei + cos 0 62 , (35a) 

= — cos 0 ei + sin (p 62 , = — sin 0 ei — cos p 62 (35b) 

with the angle 0 = 7 r /3 between slip planes and the film surfaces. The distance between crystallo¬ 
graphic slip planes is set to As = 0.05,0.1,0.2 pm, respectively. For the thickness of the crystallo¬ 
graphic layer Bs,g we choose h = 0.5As. In Case 2, we set As = 0.2, 0.4 pm. As material we use 
aluminum with a Young’s modulus oi E = 1 ■ 10^° Pa, Poisson ratio u = 0.3, Burgers vector size 
b = 2.56 • 10“^° m and drag coefficient B = 2.0 • 10“^ Pa-s. 

Numerical aspects. We use two different finite element meshes: the mesh for the elastic problem 
consists of triangular linear finite elements with altogether ~ 223 000 degrees of freedom (dofs), the 
hdCDD mesh consists of linear Fourier elements with, e.g. ~ 262 000 dofs for Case 1 with As = 100 nm 
and h = 50 nm. For the time integration we used a step size of At = 10“^ps with M = 10 ’micro-time 
steps’ per macroscopic displacement increment, resulting in 6 • 10 ^ micro-time steps for reaching the 
total strain of = 1.2%. For the shear test the same step size is used with a total number of 5 • 10^ 
micro-time steps for obtaining the total strain of = 2.5%. During each of the micro-time steps 
the stresses from the elastic BVP are kept constant and only the dislocation microstructure with the 
respective short-range interaction stresses evolves. 


12 







































































Boundary conditions. For the two different systems we introduce different boundary conditions for 
the elasticity problem: 

(1) For the symmetric tensile test we consider prescribed boundary displacements along the left 
and right boundary face (at xi = 0 and xi = 101, respectively). We increase the displacements 
with a constant rate fti^nght = —= 1.0 m/s for t G [0,0.06] ps. In order to avoid vertical 
translations we additionally fix the displacements at the point (51, 0)^. 

(2) For the shear test we consider prescribed boundary displacements along the upper surface at X 2 = 
I and fixed displacements at the lower surface at X 2 = 0. The upper prescribed displacements 
are increased with a constant rate hi,up = TO m/s for t G [0,0.05] ps. 

For both systems, also the boundary conditions for the dislocation problem w.r.t. dislocation fluxes 
have to be considered. In physical terms surfaces can either be open (dislocations can leave the 
film) or impenetrable (dislocations can not leave the film). Open boundaries can simply be modeled 
by extrapolating the hdCDD field values. For impenetrable surfaces we require (i) that the flux of 
dislocation density normal to the surface vanishes and that (ii) dislocations directly at the surface 
must be straight and thus must have zero curvature. Numerically, we model the impenetrable flux 
boundary condition by introducing a numerical inflow defined as the negative outflow of density and 
curvature densities on the considered boundary. 

Initial values. We construct consistent initial values which guarantee that, e.g., the solenoidality of 
(i.e. diva^ = 0) is not violated and that the GND density vector comes out as a gradient of 
the plastic slip. This is done by superposition of randomly distributed discrete dislocation loops 
in a 2D slip plane followed by an appropriate ’smearing-out’ procedure as described in detail in the 
appendix. One-dimensional slip plane data are then obtained by integrating the ODD field values 
over the second, homogeneous direction (see Fig. |^. Depending on the used averaging, i.e. Case 1 or 
2, we finally have to consider the distance of the representative slip planes. 


y 

Figure 3: Visualization of the initial dislocation density in the higher-dimensional configuration space 
for a system with 25 representative slip planes, each of them containing 28 ’smeared-out’ randomly 
positioned loops. In vertical direction the variable ip G [0,27r] is displayed for some of the 25 slip 
planes; the dashed white lines indicate the line orientation for screw and edge segments (compare 
Burgers vector direction in the spatial plane). Each of the wavy distribntions corresponds to (a) 
dislocation loop(s). On the bottom plane the spatial projection of the total density is shown as 
spatial average as seen in ’Case 2’. 
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4 ..!. Study 1: The influence of boundary conditions (Case 1) 

For this investigation we use 80 representative slip planes and the shear strain extension of Case 
1, starting in each SP with 5 dislocation loops of radius r between 100 nm and 200 nm at random 
positions. The height of the quasi-discrete numerical SPs has a relatively small value oi h = 50 nm 
below which no appreciable difference in the system response could be observed (also see Study 2). 
For averaging purposes we assume an out-of-plane length of = 1/ sin(0) = 1.15 pm resulting in 
an average dislocation density (p) = 3.1 x 10^^/m^. For the analysis of the influence of boundary 
conditions we study two conflgurations: in the first configuration we choose open boundaries (abbre¬ 
viated as ’open PCs’), i.e. dislocations can leave the volume, and the second imposes impenetrable 
boundaries (abbreviated as ’imp. PCs’), i.e. dislocations can not leave the volume through the surface 
dB. The simulation is driven by a prescribed constant strain rate i. = 0.2 ps“^ until the maximum 
total strain = 1.2% is reached. In Fig. |^the evolution of the total density at three distinct 


. ii 

open boundary condition ' I 

0 


2.45 


10^^/m^ impenetrable boundary condition 







\ ., y 



% A 'V ' % 




Figure 4: Initial distribution and time evolution of total density Case 1, with open and impene¬ 
trable boundary conditions for the dislocation density. 


time steps in the x — y-plane is illustrated. 

We observe that the configuration with open boundaries approaches a constant total density distri¬ 
bution along the slip planes while the system with impenetrable boundaries forms pile-ups of disloca¬ 
tions at the boundary with constant density values in between. To analyze this behavior we investigate 
the dislocation microstructure in the higher-dimensional configuration space for a slip plane in the 
center of the film in more detail (Fig. [^. The higher-dimensional fields show that for both dislocation 
boundary conditions after an initial ’incubation’ time the center region of the slip plane approaches 
a state which is characterized by only screw dislocations of positive and negative orientation: p{x^ ip) 
is approximately non-zero only for the orientations (p = 7 r /2 and p = 3 tt/2. The reason for this is 
that dislocation loops expanded and segments with edge orientation either left the film through the 
surface or pile up against the surface. In any case, only screw segments are left behind which in this 
2D model thread the film into the out-of-plane direction. Investigating the curvature k = p/q we also 
see that the screw segments are nearly straight (i.e. fc ~ 0); only for the impenetrable boundaries 
we find a non-zero curvature shortly before the surface: here, dislocations need to bend strongly in 
order to adjust from the threading screw dislocation orientation to the geometry of the Alms’ surface. 
This also suggests that the amount of dislocations inside the film for the impenetrable system will be 
significantly higher: to begin with, dislocations are not ’lost’ by out-flux through the surfaces, and an 
additionally increased line length production will take place due to the high dislocation curvature near 
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Figure 5: Time evolution of hdCDD density p(^, (y?), curvature density (p) and curvature A:(^, ip) = 
/pi^iP) for open (left block) and impenetrable boundaries (right block) for a slip plane with 
random initial values taken from the center region of Fig. On the vertical axis is the line orientation 
ip G [0, 27r], on the horizontal axis is the local ^ coordinate, the small text label indicates the maximum 
field value. 


the surfaces. Plotting the average density evolution in Fig. (b) shows that in the elastic regime I 
the dislocation density is constant, i.e. the resolved shear stress is not large enough to overcome the 
yield stress. This is followed in regime II by a transition of ’free loop expansion’ which results in 
a high dislocation multiplication rate. Towards regime III, edge components are then lost through 
open surfaces, while for impenetrable BCs edge dislocations are deposited at the surfaces. The open 
system contains at final strain a density which is smaller roughly by a factor of 3, and additionally the 
average density even reaches a stationary state (threading screw segments are straight and thus only 
translate). For the system with impenetrable surfaces the average density increases approximately 
linearly (caused by the constant line length increase of deposited edges). 

What are the consequences for the macroscopic stress-strain response? A higher dislocation density, 
on the one hand, obviously comes with a stronger influence of the Taylor equation for the yield stress. 
On the other hand, a higher plastic activity, where the density comes in through the Orowan relation, 
causes plastic softening through the solution of the elastic eigenstrain BVP. The competition of these 
two effects can be observed in the macroscopic stress-strain curve in Fig. (a) where at larger strains 
the obtained stress level for the system with open boundaries is only slightly lower. Interesting to note 
is also the initial ’hump’ right after the elastic regime I when softening sets in. This is caused by the 
fact that at an early stage dislocation loops are still comparatively small and thus contribute with a 
high line tension effect which effectively reduces the resolved shear stress. Once loops have expanded, 
this contribution is reduced and the softening behavior is sustained. A very similar behavior is also 
observed in DDD simulations (Weygand and Gumbsch, 2005[ ). 

Finally, we compare the evolution of the initial dislocation loop distribution to a distribution of 
positive and negative edge dislocations (dashed lines in Fig. |^, which would be the equivalent of a 
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Figure 6: Macroscopic quantities (Case 1) for two different boundary conditions (open and impene¬ 
trable) and two different initial configurations (consisting of either only loops or only straight edge 
dislocations). The three different time snapshots in Fig. [^correspond to the different regimes I, II, 
and III, respectively. 


’Groma-type model’. Therefore, we choose the number of edge dislocations such that the initial density 
is similar to the saturation density of loops with open BCs. We observe that after the onset of plastic 
yield edge dislocations flow towards the surfaces (regime II). They either get lost through the surfaces 
(open BCs) or pile up against the surface (imp. BCs) in which case the density simply stays constant 
because the number of edge dislocations in preserved. This microstructure results in a dramatically 
different stress-strain response as compared to the system with curved dislocations: because we have 
neither an increase in density nor a line tension one can only observe a linear hardening (regime II) 
which - regardless the boundary condition - is followed by a nearly elastic regime III where the 
stresses are considerably different from those reached for the system with dislocation loops. The loss 
of dislocations through surfaces (and ultimately the ’dislocation starvation’) has been experimentally 


observed in nano/micro pillar compression tests (Shan et ah, 2008; Jerusalem et ah, 2012) as well. 


This starvation effect also gave rise to a second elastic regime as observed in our edge dislocation 
system. In a 3D geometry we also would find this for a distribution of dislocation loops and open 
BCs. 


4-2. Study 2: Spatial coarsening from Case 1 to Case 2 

In Case 1 we considered slip planes as quasi-discrete objects mimicking the situation in DDD 
simulations. This not only requires a very high spatial resolution for the finite element scheme, but 
it is also somewhat unsatisfying from a conceptnal point of view to have two different resolutions 
(within the slip plane and perpendicular to it). If this can be avoided by use of Case 2 and whether 
it is admissible will be studied subsequently: for a given initial distribution of dislocations loops we 
compare the asymptotic system response for h ^ 0 in Case 1 and then compare with results obtained 
for different values of As in Case 2. 

All system and geometry properties are the same as in the previous Study 1, but this time we only 
compute stresses for a given dislocation configuration. To make the configurations easier to compare 
we simply choose a homogeneous distribution of loops. Fig. (a) shows the resulting stresses for the 
quasi-discrete Case 1 for a SP height of h = 50 nm and two different discretizations for Case 2. If we 
average these stress distributions we obtain stress profiles as shown in Fig. (b). It shows that for 
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(a) Distribution of the stress component <Txx 
on the left half part (0, 51) x (0, 1). 



Figure 7: Comparison of Case 1 and Case 2 for different numerical SP heights at t = 0 for a 
homogeneous distribution of loops of the same radius, (a) shows the spatial stress distribution, (b) 
shows averaged stress profiles across the height. Since the conhguration is symmetric, only half of the 
profile is shown. 


values of h = 100 nm and below the stresses do not change appreciably anymore. Running the same 
simulations for Case 2 where the plastic strain is coarse grained in between the numerical SPs and 
comparing with the results for Case 1 allows for the conclusion that for this system a As = 200 nm 
is sufficient; also the stress-strain behavior (not shown) does not show any significant difference. The 
considered situation of a homogeneous loop distribution is of course artificial. In fact, differences 
during time evolution in particular between the very coarse Case 2 and the fine Case 1 become larger 
if we start with the random initial values and as the plastic slip becomes more heterogeneous (early 
in regime II). Nonetheless, even there, our chosen approximation of Case 2 with As = 200 nm is 
sufficient. 

The advantage of the interpolation approach used in Case 2 becomes obvious when we take a look 
at the degrees of freedoms and the computational time used for the simulations from Study 1 shown 
in Tab. using Case 2 with As = 200 nm instead of Case 1 even with As = 100 nm, h = 50 nm gives 
already a speedup factor of nearly 2, while at the same time no appreciable differences in the results 
- also for the time dependent simulations - could be observed. 


configuration 

dofs (FEM + DG) 

comp, time in h on 8 procs. 

Case 1, As = 50 nm, h = 25 nm 

214962 + 1049 600 

7:39:33 

Case 1, As = 100 nm, h = 50 nm 

54682 + 262 400 

1:49:01 

Case 1, As = 200 nm, h = 100 nm 

54682 + 131200 

1:00:35 

Case 2, As = 200 nm 

54682 + 131200 

1:00:44 

Case 2, As = 400 nm 

54682 + 65 600 

0:46:04 


Table 1: Comparison of computational time and degrees of freedom for Case 1 and Case 2. 
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4-3. Study 3: A double slip configuration with Case 2 

We now investigate the macroscopic elasto-plastic response together with the microstructural evo¬ 
lution in a configuration with the two slip systems (35). We use fully averaged dislocation distributions 
(Case 2 with As = 200 nm), and the interaction of the dislocation densities in the two systems is de¬ 
scribed by the yield stress 


s,g 

^tot 


= pY' + p'r I 


(36) 


where is reconstructed in B from by averaging and interpolation using the construction 
This particular form of the yield stress represents that dislocations on the one slip system act as forest 


for the family of dislocations on the other slip system and vice versa (see e.g. Kubin et al., 2008). The 


line tension is obtained in full analogy to Study 1 for each slip system separately. Subsequently, we 
compare a single slip and double slip scenario with the following initial distribution of the dislocation 
density: in both cases we have altogether 800 dislocation loops in an averaging volume for which we 
again assume an out-of-plane averaging length of = 1.15 pm. The loops’ radii are taken from a 
uniform random distribution in the range of [100,200] nm. We either distribute them across the two 
slip systems or only across a single slip system such that the average dislocation density in the full 
body is always the same. In this study we only consider open boundaries since we focus here on the 
investigation of the hardening/softening effects introduced by the interaction of the slip systems and 
the concomitant change in dislocation microstructure. The results are illustrated in Fig. |8|-[T0 


Initially (£*°* < 0.2%) both systems respond nearly perfectly elastic because the resolved shear 
stress is in almost all regions smaller than the yield stress. This can be seen in the linear increase in 
Fig. 9a which is a consequence of the (nearly) zero velocity in Fig. (c) and (f). 

Eventually (0.2% < < 0.5%), starting from the outer regions of the density distribution 

where is smaller, the yield stress will be overcome. This results in a non-zero dislocation velocity 
mainly in the surface-near regions (Fig. |^(c) and (f) and further plots shown in appendix B). Already 
shortly after = 0.2% it becomes visible that the single slip and double slip systems behave very 
differently. The reason for this lies in the crystallography of the model systems: the Burgers vectors 
of the symmetrically inclined slip systems are such that under the prescribed shear deformation the 
diagonal components in the plastic strain tensors will cancel out. This results in a higher resolved 
shear stress since the plastic softening contribution is smaller. At the same time, however, the resulting 
velocity is higher, giving rise to more dislocation activity which can be observed in the plastic strain 
profile (compare Fig. (b) and (e)). For the single slip situation these relations are just the other 
way around: the plastic slip reduces the resolved shear stress and thus the dislocation velocity. The 
reduced plastic activity also shows in the evolution of the dislocation density and plastic slip which 
happens at a lower rate than for the double slip system, cf. Fig. 

We will now take a closer look at details of the stress-strain curve (Fig. [9a] ). What causes the 
’humps’ and different maxima for single/double slip following the elastic regime at ~ 0.5%? 


Fig. shows the higher-dimensional density, curvature density and curvature fields. There it can 
be seen, that for the single slip system the lower velocity broadens the density distribution but does 
not allow for a more significant expansion of loops and thus retards the density production. As a 
consequence of the reduced loop expansion, the loops’ curvature is also much higher for the single slip 
system, giving rise to a more pronounced infl uence of the line tension (cf. plots in Appendix B); as a 
consequence of the retarded density production the yield stress is effectively lower. The latter shows in 
Fig. [^ in the lower maximum of the hump as compared to the double slip situation; the former shows 
in the steeper inclination following the hump. While loops expand the influence of the line tension 
becomes smaller and tends to zero with edge segments leaving the film and screw dislocations becoming 
straight lines. The influence of the yield stress approaches a constant value once the density from the 
nearly straight edge dislocations saturates. Therefore, towards larger total strains (> 2%) the only 
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profiles for the single slip configuration 



(a) total density 



vertical position [nm] 


(b) plastic slip 



vertical position [nm] 


(c) velocity 


profiles for the double slip configuration 



(d) total density 



vertical position [nm] 


(e) plastic slip 



(f) velocity 


Figure 8: Profiles of CDD field variables along a central slip plane for the single slip (a-c) and double 
slip configuration (d-f) for 4 different time steps. Double slip density values and plastic slip (d-f) were 
multiplied by a factor of 2 to make them comparable with the single slip situation. 


difference between the single and double slip system stems from the plastic strain, while microscopic, 
short-range stresses had a stronger influence on early details of the stress-strain response. 

5. Summary and Conclusion 

Modeling and prediction of crystal-plasticity on the micro-meter scale requires a faithful represen¬ 
tation of the underlying physical mechanisms. We introduced the higher-dimensional CDD theory as 
a mathematical description of the kinematic behavior of statistically averaged ensembles of disloca¬ 
tions. A special emphasis was put on a consistent geometric description of the CDD field equations 
for the general case of arbitrarily oriented slip systems. In particular the transfer of information from 
two-dimensional crystallographic slip planes to the three-dimensional continuous body was concisely 
formulated. Furthermore, a conservative discontinuous Galerkin scheme suitable for the dislocation 
problem was derived. These model formulations were then applied to simulate different plain-strain 
slip geometries under tensile and shear loading conditions together with different physical boundary 
conditions. We analyzed in detail how systems of dislocation evolve while they interact with each 
other due to short-range and long-range stress components. In doing so we could directly link the 
dislocation microstructnre - represented by the higher-dimensional density and curvature density - to 
the macroscopic behavior (e.g. stress-strain response or average density). By comparing to systems of 
straight edge dislocations we observed that there the line curvature has a strong influence on the initial 
hardening behavior which is then complemented by a sustained influence of mobile screw dislocations: 
the quantitative and qualitative change in the stress-strain curve could be well explained by hdCDD. 
Furthermore, we showed how physical boundary conditions influence the orientation distribution of 
dislocations and again, how this impacts the system response. Finally, we studied a double slip con¬ 
figuration where we could attribute the distinctly different hardening behavior for single/double slip 
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(a) average resolved shear stress vs. total strain 


-- 1 slip system 

- - (7), 1 slip system 

— 2 slip systems 

— (7), 2 slip system 



total strain £*°* [%] 

(b) average total density vs. total strain 


Figure 9: Stress-strain plot and evolution of total density and plastic slip for the single slip and 
donble slip configuration. 


to microstructural aspects. 

The hdCDD model contains microstructural information which otherwise is only available in DDD 
simulations. These microstructural details have a strong influence on the system response, which is 
why a comparison with other models is difficult: no other standard continuum plasticity model is 
able to represent e.g. the conversion of SSDs into GNDs and the dislocation line length production 
accompanied by expansion of dislocation loops. For this, direct comparisons with DDD simulation 
will be extremely interesting and helpful to further benchmark our model. The necessary extraction of 
information from DDD simulations and conversion into continuous fields, however, is a non trivial task 
which needs to be nndertaken with care in future work. More realistic CDD systems will additionally 
require to incorporate further dislocation interactions and reactions as well as mechanisms for dislo¬ 


very useful there as well. With this we hope to go beyond what up to date is possible with DDD and 
apply CDD to realistic three-dimensional systems with large numbers of interacting dislocations and 
plastic strain. 
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Appendix A: Computation of consistent initial values 

Special care has to be taken with creating initial values for CDD simulations. ’Consistent initial 
values’ are those which represent averaged systems of dislocations snch that the solenoidality of a® 
is not violated, i.e. diva^ = 0. An additional constraint is the compatibility between a resulting 
GND density and the plastic shear, k-*- = — One way of guaranteeing these conditions is to 

create initial values in a 2-dimensional slip plane by superposition of objects which a priori fulfill these 
conditions. E.g. one may choose a closed dislocation loop of radius R with center point Fq; a point 


cation sources or annihilation which are not a priori included in the CDD theory (see e.g. Sandfeld 


and Hochrainer (2011) and Zhu et al. (2014) for steps into this direction). DDD simulations will be 
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Figure 10: Time evolution of hdCDD density p(^, (/?), curvature density q{^, ip) and curvature k{^, ip) = 
qi^,'p)/pi^,'p) for open boundaries, left: in single slip configuration (compare Study 1, Fig. and 
right: in double slip configuration with random initial values taken from the center region of Fig. 

On the vertical axis is the line orientation ip G [0, 27r], on the horizontal axis is the local ^ coordinate, 
the small text label indicate the maximum field value (comparison with Fig. [^of Study 1). 


r in Fs^g then has the distance d = |i?(—sine/?,cos(y?) + Fq — r| from the loops’ line and the density 
and the curvature can be obtained by ’smearing-out’ the line such that the total line length 27rR is 
preserved, i.e., 


p(r,<F) 


fiexp(-l/(l-(d/do)^)) 

2n texp{-l/{l-{t/doP))dt 

0 


if d < do , 
else. 


where the parameter do governs the width of the compact density distribution around the line. Note, 
that this resnlts in an area density i.e. p has the unit of line length per area. The curvature density 

can be obtained as q{v,p) = The corresponding plastic slip generated by the expanded 

pR-\-do p27t 

loop is given by 7 (r) =b / p{r], 0, p) dp dp. We superimpose such dislocation loops with 


^|r-ro| -'0 

randomly chosen centers of the loops with a uniform distribution and such that the compact function 
for the CDD fields fits completely into the slip plane. Additionally, the loops’ radii are chosen from 
a uniform random distribution. The ID-slip planes represent a homogeneous distribution into the 
second direction, so that we integrate the CDD fields along the width for Ps,gi0^Qs,g{0 
see Fig. |^for an visualization. 


Appendix B: Evolution of CDD field quantities for Study 3 

The following figures show the evolution of total density for the single and double slip system 
(Fig. ng. The largest difference in these two situations is that due to the different stress state in 
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double slip the dislocation activity is much more pronounced. In Fig.[T^also the profiles for all relevant 
CDD field quantities and stress components are shown. One of the key feature of hdCDD is that the 
distinction into GNDs and SSDs is obtained naturally from the higher-dimensional configuration space. 
This can be also observed in the evolution of the integrated quantities and in Fig. 
and (g)+(h). Therein, the difference would yield the SSD density of edges. 
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Figure 11: Initial distribution and time evolution of total density pl^g for Case 2, Study 3, for system 
with 5 = 1 and s = 1,2 at t = 0.01 ps and t = 0.02 ps with open boundary conditions for the 
dislocation density. 
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